##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("effectsize")) install.packages("effectsize")
if (!require("emmeans")) install.packages("emmeans")
Loading required package: emmeans
if (!require("dplyr")) install.packages("dplyr")
if (!require("tidyr")) install.packages("tidyr")
Loading required package: tidyr
if (!require("Rmisc")) install.packages("Rmisc")
Loading required package: Rmisc
Loading required package: plyr
----------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------------------------------------
Attaching package: ‘plyr’
The following objects are masked from ‘package:dplyr’:
arrange, count, desc, failwith, id, mutate, rename, summarise,
summarize
if (!require("ggpubr")) install.packages("ggpubr")
Loading required package: ggpubr
Attaching package: ‘ggpubr’
The following object is masked from ‘package:plyr’:
mutate
The following object is masked from ‘package:cowplot’:
get_legend
if (!require("cowplot")) install.packages("cowplot")
if (!require("gridGraphics")) install.packages("gridGraphics")
Loading required package: gridGraphics
Warning: package ‘gridGraphics’ was built under R version 4.3.2Loading required package: grid
if (!require("tidyr")) install.packages("tidyr")
##Load Packages
library(ggplot2) #Required for plotting
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons
library(dplyr) #Required to seperate columns in dataframe
library(tidyr) #Required for data organization
library(Rmisc) #Required for summarySE for summary statistics
library(ggpubr) #Required for adding pairwise p-values to plots with stat_pvalue_manual
library(cowplot) #Required for arranging ggplots
library(gridGraphics) #Required for adding labels to arranged plots
library(tidyr) #Required for reshaping datafrom from a wide to long format.
Note: Run “Graphing Parameters” section from 01_ExperimentalSetup.R file
Note: Full Data with Bleaching Metrics and Color Scores created in 04_Models.R file
#Load Data
FullData<-read.csv("Outputs/FullData.csv", header=TRUE)
#Set factor variables
FullData$TimeP<-factor(FullData$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
FullData$Site<-factor(FullData$Site, levels=c("KL", "SS"), ordered=TRUE)
FullData$Genotype<-factor(FullData$Genotype, levels=c("AC10", "AC12", "AC8"), ordered=TRUE)
FullData$Treatment<-factor(FullData$Treatment, levels=c("Control", "Heat"), ordered=TRUE)
FullData$Treat<-factor(FullData$Treat, levels=c("C", "H"), ordered=TRUE)
FullData$Seas<-factor(FullData$Seas, levels=c("Summer1", "Summer2", "Winter"), ordered=TRUE)
##Control
FullData_C<-subset(FullData, Treat=="C")
##Heated
FullData_H<-subset(FullData, Treat=="H")
Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.
##Calculate averages for Control Treatment for each Site, Genotype, Timepoint, and Season
names(FullData_C)
[1] "ID" "RandN" "TimeP" "Site" "Genotype"
[6] "Treat" "Treatment" "Seas" "Set" "Score_Full"
[11] "Score_TP" "SA_cm2" "Chl_ug.cm2" "Sym10.6_cm2"
FullData_C.a<-aggregate(FullData_C[,c(10:11, 13:14)], list(FullData_C$Site, FullData_C$Genotype, FullData_C$TimeP, FullData_C$Seas), mean, na.action = na.omit)
names(FullData_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Seas")
names(FullData_C.a)[5:8]<-paste(names(FullData_C)[c(10:11, 13:14)], "C", sep="_")
##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, Timepoint, and Season
TolData<-merge(FullData_H, FullData_C.a, all.x=TRUE )
##Calculate Proportion Retained for each Bleaching Metric relative to Control
TolData$Score_Full.prop<-round((TolData$Score_Full/TolData$Score_Full_C), 4)
TolData$Score_TP.prop<-round((TolData$Score_TP/TolData$Score_TP_C), 4)
TolData$Chl.prop<-round((TolData$Chl_ug.cm2/TolData$Chl_ug.cm2_C), 4)
TolData$Sym.prop<-round((TolData$Sym10.6_cm2/TolData$Sym10.6_cm2_C), 4)
##Set values >1 to 1
TolData$Score_Full.prop[which(TolData$Score_Full.prop>1)]<-1.0000
TolData$Score_TP.prop[which(TolData$Score_TP.prop>1)]<-1.0000
TolData$Chl.prop[which(TolData$Chl.prop>1)]<-1.0000
TolData$Sym.prop[which(TolData$Sym.prop>1)]<-1.0000
##Create Site and Genotype Variable
TolData$Site.Geno<-paste(TolData$Site, TolData$Genotype, sep="_")
TolData_W2<-subset(TolData, TimeP=="W2")
TolData_M1<-subset(TolData, TimeP=="M1")
TolData_M4<-subset(TolData, TimeP=="M4")
##Check normality
hist(TolData$Chl.prop)
shapiro.test(TolData$Chl.prop)
Shapiro-Wilk normality test
data: TolData$Chl.prop
W = 0.97395, p-value = 0.165
#Normal
##Model as a function of Site and Genotype and Timepoint
Tol_Chl_lm<-lm(Chl.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_lm)); qqline(resid(Tol_Chl_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_lm), resid(Tol_Chl_lm))
#Model Results
summary(Tol_Chl_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.127060 -0.031006 -0.004219 0.036829 0.097690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2100753 0.0065828 31.913 < 2e-16 ***
Site.L -0.0093229 0.0092942 -1.003 0.320289
Genotype.L -0.0980285 0.0111500 -8.792 5.31e-12 ***
Genotype.Q -0.0467859 0.0116301 -4.023 0.000180 ***
TimeP.L 0.0638187 0.0112994 5.648 6.22e-07 ***
TimeP.Q 0.0464233 0.0115033 4.036 0.000173 ***
Site.L:Genotype.L 0.0003917 0.0157685 0.025 0.980272
Site.L:Genotype.Q 0.0592529 0.0164183 3.609 0.000673 ***
Site.L:TimeP.L -0.0528050 0.0159798 -3.304 0.001694 **
Site.L:TimeP.Q -0.0236749 0.0161894 -1.462 0.149435
Genotype.L:TimeP.L 0.0406686 0.0194339 2.093 0.041093 *
Genotype.Q:TimeP.L -0.0627603 0.0196595 -3.192 0.002354 **
Genotype.L:TimeP.Q 0.0800244 0.0191901 4.170 0.000111 ***
Genotype.Q:TimeP.Q -0.1328570 0.0206170 -6.444 3.28e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05393 on 54 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.8229, Adjusted R-squared: 0.7803
F-statistic: 19.3 on 13 and 54 DF, p-value: 8.612e-16
anova(Tol_Chl_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.001563 0.001563 0.5375 0.466633
Genotype 2 0.297307 0.148654 51.1117 3.495e-13 ***
TimeP 2 0.132285 0.066142 22.7418 6.845e-08 ***
Site:Genotype 2 0.041114 0.020557 7.0681 0.001877 **
Site:TimeP 2 0.038643 0.019322 6.6434 0.002634 **
Genotype:TimeP 4 0.218951 0.054738 18.8205 9.691e-10 ***
Residuals 54 0.157054 0.002908
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 1.76e-03 | [0.00, 1.00]
Genotype | 0.34 | [0.16, 1.00]
TimeP | 0.15 | [0.02, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
Site:TimeP | 0.04 | [0.00, 1.00]
Genotype:TimeP | 0.25 | [0.06, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_lm.res<-data.frame(anova(Tol_Chl_lm))
Tol_Chl_lm.res$Predictor<-rownames(Tol_Chl_lm.res)
Tol_Chl_lm.res$EtaSq<-c(eta_squared(Tol_Chl_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_lm.res))
Tol_Chl_lm.res<-Tol_Chl_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Chl.prop)
shapiro.test(TolData_W2$Chl.prop)
Shapiro-Wilk normality test
data: TolData_W2$Chl.prop
W = 0.98407, p-value = 0.9671
#Normal
##Model as a function of Site and Genotype
Tol_Chl_W2_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_W2_lm)); qqline(resid(Tol_Chl_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_W2_lm), resid(Tol_Chl_W2_lm))
#Model Results
summary(Tol_Chl_W2_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.089900 -0.031587 -0.008183 0.038881 0.103550
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1842264 0.0119450 15.423 5.03e-11 ***
Site.L 0.0188110 0.0168928 1.114 0.281921
Genotype.L -0.0947759 0.0204291 -4.639 0.000273 ***
Genotype.Q -0.0585870 0.0209464 -2.797 0.012921 *
Site.L:Genotype.L -0.0002833 0.0288911 -0.010 0.992297
Site.L:Genotype.Q 0.0366569 0.0296227 1.237 0.233769
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05552 on 16 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.6758, Adjusted R-squared: 0.5745
F-statistic: 6.671 on 5 and 16 DF, p-value: 0.001551
anova(Tol_Chl_W2_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.005580 0.005580 1.8106 0.197195
Genotype 2 0.092482 0.046241 15.0038 0.000214 ***
Site:Genotype 2 0.004732 0.002366 0.7677 0.480448
Residuals 16 0.049311 0.003082
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.04 | [0.00, 1.00]
Genotype | 0.61 | [0.29, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_W2_lm.res<-data.frame(anova(Tol_Chl_W2_lm))
Tol_Chl_W2_lm.res$Predictor<-rownames(Tol_Chl_W2_lm.res)
Tol_Chl_W2_lm.res$EtaSq<-c(eta_squared(Tol_Chl_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_W2_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res<-Tol_Chl_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.2033 0.0278 16 0.1445 0.262
AC12 0.2399 0.0278 16 0.1811 0.299
AC8 0.0696 0.0278 16 0.0107 0.128
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.2514 0.0278 16 0.1925 0.310
AC12 0.2242 0.0321 16 0.1563 0.292
AC8 0.1170 0.0321 16 0.0491 0.185
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.0366 0.0393 16 -0.933 0.6280
AC10 - AC8 0.1338 0.0393 16 3.407 0.0095
AC12 - AC8 0.1704 0.0393 16 4.340 0.0014
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0272 0.0424 16 0.640 0.8003
AC10 - AC8 0.1343 0.0424 16 3.168 0.0156
AC12 - AC8 0.1072 0.0453 16 2.364 0.0753
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.2033 0.0278 16 0.1445 0.262
SS 0.2514 0.0278 16 0.1925 0.310
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.2399 0.0278 16 0.1811 0.299
SS 0.2242 0.0321 16 0.1563 0.292
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0696 0.0278 16 0.0107 0.128
SS 0.1170 0.0321 16 0.0491 0.185
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0481 0.0393 16 -1.224 0.2387
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.0157 0.0424 16 0.371 0.7156
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0475 0.0424 16 -1.120 0.2793
##Save p-values
#Genotypes within Sites
Tol_Chl_W2_lm.geno<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_W2_lm.geno<-Tol_Chl_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.geno$group1<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group1, sep="_")
Tol_Chl_W2_lm.geno$group2<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_W2_lm.site<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_W2_lm.site<-Tol_Chl_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.site$group1<-paste(Tol_Chl_W2_lm.site$group1, Tol_Chl_W2_lm.site$Genotype, sep="_")
Tol_Chl_W2_lm.site$group2<-paste(Tol_Chl_W2_lm.site$group2, Tol_Chl_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_W2_lm.p<-rbind(Tol_Chl_W2_lm.geno[,c(1:2,4:8)], Tol_Chl_W2_lm.site[,c(1:2,4:8)])
Tol_Chl_W2_lm.p<-Tol_Chl_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_W2_lm.p$Sig<-ifelse(Tol_Chl_W2_lm.p$p<0.001, "***", ifelse(Tol_Chl_W2_lm.p$p<0.01, "**", ifelse(Tol_Chl_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_W2_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.p))
Tol_Chl_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_W2_SG<-summarySE(TolData_W2, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_W2_SG.plot<-ggplot(Tol_Chl_W2_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Chlorophyll Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_W2_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_W2_SG.plot
##Check normality
hist(TolData_M1$Chl.prop)
shapiro.test(TolData_M1$Chl.prop)
Shapiro-Wilk normality test
data: TolData_M1$Chl.prop
W = 0.9306, p-value = 0.1263
#Normal
##Model as a function of Site and Genotype
Tol_Chl_M1_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_M1_lm)); qqline(resid(Tol_Chl_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M1_lm), resid(Tol_Chl_M1_lm))
#Model Results
summary(Tol_Chl_M1_lm)
Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.107050 -0.034575 0.003479 0.025531 0.088050
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.17217 0.01178 14.619 1.12e-10 ***
Site.L 0.01036 0.01666 0.622 0.5428
Genotype.L -0.16337 0.01935 -8.442 2.74e-07 ***
Genotype.Q 0.06169 0.02139 2.883 0.0108 *
Site.L:Genotype.L 0.03819 0.02737 1.395 0.1820
Site.L:Genotype.Q 0.05454 0.03026 1.803 0.0903 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05474 on 16 degrees of freedom
Multiple R-squared: 0.8424, Adjusted R-squared: 0.7932
F-statistic: 17.11 on 5 and 16 DF, p-value: 6.411e-06
anova(Tol_Chl_M1_lm)
Analysis of Variance Table
Response: Chl.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.002283 0.002283 0.7619 0.3956
Genotype 2 0.238424 0.119212 39.7888 6.168e-07 ***
Site:Genotype 2 0.015569 0.007785 2.5982 0.1054
Residuals 16 0.047938 0.002996
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 7.50e-03 | [0.00, 1.00]
Genotype | 0.78 | [0.58, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_M1_lm.res<-data.frame(anova(Tol_Chl_M1_lm))
Tol_Chl_M1_lm.res$Predictor<-rownames(Tol_Chl_M1_lm.res)
Tol_Chl_M1_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M1_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res<-Tol_Chl_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.3089 0.0274 16 0.2509 0.3669
AC12 0.1460 0.0316 16 0.0790 0.2130
AC8 0.0397 0.0274 16 -0.0183 0.0977
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.3169 0.0274 16 0.2588 0.3749
AC12 0.0976 0.0316 16 0.0306 0.1646
AC8 0.1240 0.0274 16 0.0660 0.1820
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1629 0.0418 16 3.897 0.0035
AC10 - AC8 0.2692 0.0387 16 6.956 <.0001
AC12 - AC8 0.1063 0.0418 16 2.543 0.0538
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.2192 0.0418 16 5.244 0.0002
AC10 - AC8 0.1928 0.0387 16 4.983 0.0004
AC12 - AC8 -0.0264 0.0418 16 -0.631 0.8056
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.3089 0.0274 16 0.2509 0.3669
SS 0.3169 0.0274 16 0.2588 0.3749
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.1460 0.0316 16 0.0790 0.2130
SS 0.0976 0.0316 16 0.0306 0.1646
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0397 0.0274 16 -0.0183 0.0977
SS 0.1240 0.0274 16 0.0660 0.1820
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00795 0.0387 16 -0.205 0.8398
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.04833 0.0447 16 1.081 0.2955
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.08432 0.0387 16 -2.179 0.0447
##Save p-values
#Genotypes within Sites
Tol_Chl_M1_lm.geno<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M1_lm.geno<-Tol_Chl_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.geno$group1<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group1, sep="_")
Tol_Chl_M1_lm.geno$group2<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_M1_lm.site<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M1_lm.site<-Tol_Chl_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.site$group1<-paste(Tol_Chl_M1_lm.site$group1, Tol_Chl_M1_lm.site$Genotype, sep="_")
Tol_Chl_M1_lm.site$group2<-paste(Tol_Chl_M1_lm.site$group2, Tol_Chl_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_M1_lm.p<-rbind(Tol_Chl_M1_lm.geno[,c(1:2,4:8)], Tol_Chl_M1_lm.site[,c(1:2,4:8)])
Tol_Chl_M1_lm.p<-Tol_Chl_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_M1_lm.p$Sig<-ifelse(Tol_Chl_M1_lm.p$p<0.001, "***", ifelse(Tol_Chl_M1_lm.p$p<0.01, "**", ifelse(Tol_Chl_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_M1_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.p))
Tol_Chl_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_M1_SG<-summarySE(TolData_M1, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_M1_SG.plot<-ggplot(Tol_Chl_M1_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Chlorophyll Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M1_lm.p, y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_M1_SG.plot
##Check normality
hist(TolData_M4$Chl.prop)
shapiro.test(TolData_M4$Chl.prop)
Shapiro-Wilk normality test
data: TolData_M4$Chl.prop
W = 0.91299, p-value = 0.04096
#Not Normal
##Try log+1 transformation
hist(log(TolData_M4$Chl.prop+1))
shapiro.test(log(TolData_M4$Chl.prop+1))
Shapiro-Wilk normality test
data: log(TolData_M4$Chl.prop + 1)
W = 0.93102, p-value = 0.1028
#Normal
##Model as a function of Site and Genotype
##Model with log transformation
Tol_Chl_M4_lm<-lm(log(Chl.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Chl_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_Chl_M4_lm)); qqline(resid(Tol_Chl_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M4_lm), resid(Tol_Chl_M4_lm))
#Model Results
summary(Tol_Chl_M4_lm)
Call:
lm(formula = log(Chl.prop + 1) ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.062753 -0.025900 -0.004896 0.017009 0.065759
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.23863 0.00807 29.568 < 2e-16 ***
Site.L -0.04272 0.01141 -3.743 0.00149 **
Genotype.L -0.03027 0.01398 -2.166 0.04400 *
Genotype.Q -0.11052 0.01398 -7.907 2.9e-07 ***
Site.L:Genotype.L -0.03189 0.01977 -1.613 0.12415
Site.L:Genotype.Q 0.05753 0.01977 2.910 0.00934 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03954 on 18 degrees of freedom
Multiple R-squared: 0.8368, Adjusted R-squared: 0.7915
F-statistic: 18.46 on 5 and 18 DF, p-value: 1.598e-06
anova(Tol_Chl_M4_lm)
Analysis of Variance Table
Response: log(Chl.prop + 1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.021901 0.021901 14.0102 0.001489 **
Genotype 2 0.105056 0.052528 33.6031 8.379e-07 ***
Site:Genotype 2 0.017304 0.008652 5.5349 0.013381 *
Residuals 18 0.028137 0.001563
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Chl_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.13 | [0.00, 1.00]
Genotype | 0.61 | [0.32, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Chl_M4_lm.res<-data.frame(anova(Tol_Chl_M4_lm))
Tol_Chl_M4_lm.res$Predictor<-rownames(Tol_Chl_M4_lm.res)
Tol_Chl_M4_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M4_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res<-Tol_Chl_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype and Site*Genotype.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.213 0.0198 18 0.171 0.254
AC12 0.392 0.0198 18 0.351 0.434
AC8 0.202 0.0198 18 0.160 0.243
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.217 0.0198 18 0.176 0.259
AC12 0.265 0.0198 18 0.224 0.307
AC8 0.143 0.0198 18 0.101 0.184
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.1797 0.028 18 -6.428 <.0001
AC10 - AC8 0.0109 0.028 18 0.391 0.9196
AC12 - AC8 0.1906 0.028 18 6.819 <.0001
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.0482 0.028 18 -1.724 0.2237
AC10 - AC8 0.0747 0.028 18 2.672 0.0393
AC12 - AC8 0.1229 0.028 18 4.396 0.0010
Note: contrasts are still on the log(mu + 1) scale
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.213 0.0198 18 0.171 0.254
SS 0.217 0.0198 18 0.176 0.259
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.392 0.0198 18 0.351 0.434
SS 0.265 0.0198 18 0.224 0.307
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.202 0.0198 18 0.160 0.243
SS 0.143 0.0198 18 0.101 0.184
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00468 0.028 18 -0.167 0.8688
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.12684 0.028 18 4.537 0.0003
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.05909 0.028 18 2.114 0.0488
Note: contrasts are still on the log(mu + 1) scale
##Save p-values
#Genotypes within Sites
Tol_Chl_M4_lm.geno<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M4_lm.geno<-Tol_Chl_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.geno$group1<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group1, sep="_")
Tol_Chl_M4_lm.geno$group2<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Chl_M4_lm.site<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M4_lm.site<-Tol_Chl_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.site$group1<-paste(Tol_Chl_M4_lm.site$group1, Tol_Chl_M4_lm.site$Genotype, sep="_")
Tol_Chl_M4_lm.site$group2<-paste(Tol_Chl_M4_lm.site$group2, Tol_Chl_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Chl_M4_lm.p<-rbind(Tol_Chl_M4_lm.geno[,c(1:2,4:8)], Tol_Chl_M4_lm.site[,c(1:2,4:8)])
Tol_Chl_M4_lm.p<-Tol_Chl_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Chl_M4_lm.p$Sig<-ifelse(Tol_Chl_M4_lm.p$p<0.001, "***", ifelse(Tol_Chl_M4_lm.p$p<0.01, "**", ifelse(Tol_Chl_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Chl_M4_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.p))
Tol_Chl_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_Chl_M4_SG<-summarySE(TolData_M4, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Chl_M4_SG.plot<-ggplot(Tol_Chl_M4_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Chlorophyll Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Chl_M4_lm.p, y.position=0.55, step.increase=0.20, label="Sig", hide.ns=TRUE); Tol_Chl_M4_SG.plot
##Check normality
hist(TolData$Sym.prop)
shapiro.test(TolData$Sym.prop)
Shapiro-Wilk normality test
data: TolData$Sym.prop
W = 0.92235, p-value = 0.0003707
#Not Normal
##Try log+1 transformation
hist(log(TolData$Sym.prop+1))
shapiro.test(log(TolData$Sym.prop+1))
Shapiro-Wilk normality test
data: log(TolData$Sym.prop + 1)
W = 0.92164, p-value = 0.000345
#Not normal
##Try square transformation
hist((TolData$Sym.prop)^2)
shapiro.test((TolData$Sym.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Sym.prop)^2
W = 0.87147, p-value = 3.933e-06
#Not normal
##Try cubed transformation
hist((TolData$Sym.prop)^3)
shapiro.test((TolData$Sym.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Sym.prop)^3
W = 0.80697, p-value = 4.342e-08
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_Sym_lm<-lm(Sym.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_lm)); qqline(resid(Tol_Sym_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_lm), resid(Tol_Sym_lm))
#Model Results
summary(Tol_Sym_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.37464 -0.08630 -0.00192 0.08972 0.35883
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.544513 0.018576 29.313 < 2e-16 ***
Site.L -0.001834 0.026197 -0.070 0.944443
Genotype.L -0.373179 0.031737 -11.758 < 2e-16 ***
Genotype.Q 0.028101 0.032605 0.862 0.392497
TimeP.L 0.133503 0.031737 4.206 9.64e-05 ***
TimeP.Q 0.190344 0.032605 5.838 2.94e-07 ***
Site.L:Genotype.L -0.060389 0.044883 -1.345 0.183999
Site.L:Genotype.Q 0.074992 0.045930 1.633 0.108240
Site.L:TimeP.L -0.159707 0.044883 -3.558 0.000778 ***
Site.L:TimeP.Q -0.027928 0.045930 -0.608 0.545653
Genotype.L:TimeP.L 0.085371 0.055316 1.543 0.128489
Genotype.Q:TimeP.L -0.084924 0.054623 -1.555 0.125749
Genotype.L:TimeP.Q 0.065986 0.054623 1.208 0.232204
Genotype.Q:TimeP.Q -0.205029 0.058264 -3.519 0.000878 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1535 on 55 degrees of freedom
Multiple R-squared: 0.8036, Adjusted R-squared: 0.7572
F-statistic: 17.31 on 13 and 55 DF, p-value: 6.022e-15
anova(Tol_Sym_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.0007 0.00068 0.0289 0.865708
Genotype 2 3.2796 1.63981 69.5880 8.602e-16 ***
TimeP 2 1.1556 0.57780 24.5197 2.438e-08 ***
Site:Genotype 2 0.1111 0.05557 2.3581 0.104096
Site:TimeP 2 0.3138 0.15688 6.6576 0.002575 **
Genotype:TimeP 4 0.4428 0.11069 4.6974 0.002474 **
Residuals 55 1.2961 0.02356
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 1.03e-04 | [0.00, 1.00]
Genotype | 0.50 | [0.33, 1.00]
TimeP | 0.18 | [0.04, 1.00]
Site:Genotype | 0.02 | [0.00, 1.00]
Site:TimeP | 0.05 | [0.00, 1.00]
Genotype:TimeP | 0.07 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_lm.res<-data.frame(anova(Tol_Sym_lm))
Tol_Sym_lm.res$Predictor<-rownames(Tol_Sym_lm.res)
Tol_Sym_lm.res$EtaSq<-c(eta_squared(Tol_Sym_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_lm.res))
Tol_Sym_lm.res<-Tol_Sym_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Sym.prop)
shapiro.test(TolData_W2$Sym.prop)
Shapiro-Wilk normality test
data: TolData_W2$Sym.prop
W = 0.93605, p-value = 0.1477
#Normal
##Model as a function of Site and Genotype
Tol_Sym_W2_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_W2_lm)); qqline(resid(Tol_Sym_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_W2_lm), resid(Tol_Sym_W2_lm))
#Model Results
summary(Tol_Sym_W2_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.26060 -0.11145 0.00000 0.08525 0.28580
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5237625 0.0325072 16.112 9.90e-12 ***
Site.L 0.0939568 0.0459722 2.044 0.0568 .
Genotype.L -0.4152131 0.0570402 -7.279 1.29e-06 ***
Genotype.Q -0.0005205 0.0555584 -0.009 0.9926
Site.L:Genotype.L -0.2061000 0.0806671 -2.555 0.0205 *
Site.L:Genotype.Q 0.0462891 0.0785715 0.589 0.5635
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.155 on 17 degrees of freedom
Multiple R-squared: 0.7924, Adjusted R-squared: 0.7313
F-statistic: 12.98 on 5 and 17 DF, p-value: 2.647e-05
anova(Tol_Sym_W2_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.16148 0.16148 6.7208 0.01897 *
Genotype 2 1.22843 0.61421 25.5640 7.508e-06 ***
Site:Genotype 2 0.16883 0.08441 3.5133 0.05283 .
Residuals 17 0.40845 0.02403
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.08 | [0.00, 1.00]
Genotype | 0.62 | [0.33, 1.00]
Site:Genotype | 0.09 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_W2_lm.res<-data.frame(anova(Tol_Sym_W2_lm))
Tol_Sym_W2_lm.res$Predictor<-rownames(Tol_Sym_W2_lm.res)
Tol_Sym_W2_lm.res$EtaSq<-c(eta_squared(Tol_Sym_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_W2_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res<-Tol_Sym_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.634 0.0775 17 0.4708 0.798
AC12 0.484 0.0775 17 0.3210 0.648
AC8 0.253 0.0775 17 0.0897 0.417
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 1.000 0.0775 17 0.8365 1.164
AC12 0.564 0.0775 17 0.4004 0.727
AC8 0.207 0.0895 17 0.0179 0.396
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.150 0.110 17 1.367 0.3795
AC10 - AC8 0.381 0.110 17 3.477 0.0077
AC12 - AC8 0.231 0.110 17 2.110 0.1174
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.436 0.110 17 3.979 0.0026
AC10 - AC8 0.793 0.118 17 6.701 <.0001
AC12 - AC8 0.357 0.118 17 3.017 0.0201
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.634 0.0775 17 0.4708 0.798
SS 1.000 0.0775 17 0.8365 1.164
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.484 0.0775 17 0.3210 0.648
SS 0.564 0.0775 17 0.4004 0.727
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.253 0.0775 17 0.0897 0.417
SS 0.207 0.0895 17 0.0179 0.396
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.3657 0.110 17 -3.337 0.0039
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0794 0.110 17 -0.725 0.4785
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.0465 0.118 17 0.393 0.6994
##Save p-values
#Genotypes within Sites
Tol_Sym_W2_lm.geno<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_W2_lm.geno<-Tol_Sym_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.geno$group1<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group1, sep="_")
Tol_Sym_W2_lm.geno$group2<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_W2_lm.site<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_W2_lm.site<-Tol_Sym_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.site$group1<-paste(Tol_Sym_W2_lm.site$group1, Tol_Sym_W2_lm.site$Genotype, sep="_")
Tol_Sym_W2_lm.site$group2<-paste(Tol_Sym_W2_lm.site$group2, Tol_Sym_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_W2_lm.p<-rbind(Tol_Sym_W2_lm.geno[,c(1:2,4:8)], Tol_Sym_W2_lm.site[,c(1:2,4:8)])
Tol_Sym_W2_lm.p<-Tol_Sym_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_W2_lm.p$Sig<-ifelse(Tol_Sym_W2_lm.p$p<0.001, "***", ifelse(Tol_Sym_W2_lm.p$p<0.01, "**", ifelse(Tol_Sym_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_W2_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.p))
Tol_Sym_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_W2_SG<-summarySE(TolData_W2, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_W2_SG.plot<-ggplot(Tol_Sym_W2_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Symbiont Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_W2_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_W2_SG.plot
##Check normality
hist(TolData_M1$Sym.prop)
shapiro.test(TolData_M1$Sym.prop)
Shapiro-Wilk normality test
data: TolData_M1$Sym.prop
W = 0.85748, p-value = 0.004618
#Not normal
##Try log+1 transformation
hist(log(TolData_M1$Sym.prop+1))
shapiro.test(log(TolData_M1$Sym.prop+1))
Shapiro-Wilk normality test
data: log(TolData_M1$Sym.prop + 1)
W = 0.87335, p-value = 0.009041
#Not normal but improved
##Model as a function of Site and Genotype
##Model with log transformation
Tol_Sym_M1_lm<-lm(log(Sym.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_M1_lm)); qqline(resid(Tol_Sym_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M1_lm), resid(Tol_Sym_M1_lm))
#Model Results
summary(Tol_Sym_M1_lm)
Call:
lm(formula = log(Sym.prop + 1) ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.098079 -0.055799 -0.004888 0.043331 0.151113
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.307716 0.016556 18.587 2.95e-12 ***
Site.L 0.027767 0.023413 1.186 0.252954
Genotype.L -0.295494 0.027204 -10.862 8.59e-09 ***
Genotype.Q 0.129387 0.030075 4.302 0.000548 ***
Site.L:Genotype.L 0.080682 0.038472 2.097 0.052224 .
Site.L:Genotype.Q -0.005572 0.042532 -0.131 0.897412
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.07694 on 16 degrees of freedom
Multiple R-squared: 0.8989, Adjusted R-squared: 0.8673
F-statistic: 28.46 on 5 and 16 DF, p-value: 2e-07
anova(Tol_Sym_M1_lm)
Analysis of Variance Table
Response: log(Sym.prop + 1)
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.00823 0.00823 1.3902 0.2556
Genotype 2 0.80811 0.40406 68.2486 1.468e-08 ***
Site:Genotype 2 0.02614 0.01307 2.2076 0.1423
Residuals 16 0.09473 0.00592
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 8.78e-03 | [0.00, 1.00]
Genotype | 0.86 | [0.73, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_M1_lm.res<-data.frame(anova(Tol_Sym_M1_lm))
Tol_Sym_M1_lm.res$Predictor<-rownames(Tol_Sym_M1_lm.res)
Tol_Sym_M1_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M1_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res<-Tol_Sym_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.5918 0.0385 16 0.5102 0.673
AC12 0.1792 0.0444 16 0.0850 0.273
AC8 0.0932 0.0385 16 0.0117 0.175
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.5472 0.0385 16 0.4656 0.629
AC12 0.2249 0.0444 16 0.1307 0.319
AC8 0.2100 0.0385 16 0.1284 0.292
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.413 0.0588 16 7.021 <.0001
AC10 - AC8 0.499 0.0544 16 9.164 <.0001
AC12 - AC8 0.086 0.0588 16 1.463 0.3339
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.322 0.0588 16 5.483 0.0001
AC10 - AC8 0.337 0.0544 16 6.198 <.0001
AC12 - AC8 0.015 0.0588 16 0.255 0.9650
Note: contrasts are still on the log(mu + 1) scale
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.5918 0.0385 16 0.5102 0.673
SS 0.5472 0.0385 16 0.4656 0.629
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.1792 0.0444 16 0.0850 0.273
SS 0.2249 0.0444 16 0.1307 0.319
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.0932 0.0385 16 0.0117 0.175
SS 0.2100 0.0385 16 0.1284 0.292
Results are given on the log(mu + 1) (not the response) scale.
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS 0.0446 0.0544 16 0.820 0.4241
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0457 0.0628 16 -0.727 0.4775
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1167 0.0544 16 -2.146 0.0476
Note: contrasts are still on the log(mu + 1) scale
##Save p-values
#Genotypes within Sites
Tol_Sym_M1_lm.geno<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M1_lm.geno<-Tol_Sym_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.geno$group1<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group1, sep="_")
Tol_Sym_M1_lm.geno$group2<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_M1_lm.site<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M1_lm.site<-Tol_Sym_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.site$group1<-paste(Tol_Sym_M1_lm.site$group1, Tol_Sym_M1_lm.site$Genotype, sep="_")
Tol_Sym_M1_lm.site$group2<-paste(Tol_Sym_M1_lm.site$group2, Tol_Sym_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_M1_lm.p<-rbind(Tol_Sym_M1_lm.geno[,c(1:2,4:8)], Tol_Sym_M1_lm.site[,c(1:2,4:8)])
Tol_Sym_M1_lm.p<-Tol_Sym_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_M1_lm.p$Sig<-ifelse(Tol_Sym_M1_lm.p$p<0.001, "***", ifelse(Tol_Sym_M1_lm.p$p<0.01, "**", ifelse(Tol_Sym_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_M1_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.p))
Tol_Sym_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_M1_SG<-summarySE(TolData_M1, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_M1_SG.plot<-ggplot(Tol_Sym_M1_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Symbiont Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M1_lm.p, y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M1_SG.plot
##Check normality
hist(TolData_M4$Sym.prop)
shapiro.test(TolData_M4$Sym.prop)
Shapiro-Wilk normality test
data: TolData_M4$Sym.prop
W = 0.88557, p-value = 0.01078
#Not Normal
##Try square transformation
hist((TolData_M4$Sym.prop)^2)
shapiro.test((TolData_M4$Sym.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Sym.prop)^2
W = 0.87973, p-value = 0.008208
#Not Normal
##Try cubed transformation
hist((TolData_M4$Sym.prop)^3)
shapiro.test((TolData_M4$Sym.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Sym.prop)^3
W = 0.8488, p-value = 0.002077
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_Sym_M4_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_Sym_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_Sym_M4_lm)); qqline(resid(Tol_Sym_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M4_lm), resid(Tol_Sym_M4_lm))
#Model Results
summary(Tol_Sym_M4_lm)
Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.31820 -0.08804 0.00000 0.08618 0.29740
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.71662 0.03157 22.702 1.07e-14 ***
Site.L -0.12617 0.04464 -2.826 0.0112 *
Genotype.L -0.28587 0.05467 -5.229 5.67e-05 ***
Genotype.Q -0.11565 0.05467 -2.115 0.0486 *
Site.L:Genotype.L -0.09494 0.07732 -1.228 0.2353
Site.L:Genotype.Q 0.17275 0.07732 2.234 0.0384 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1546 on 18 degrees of freedom
Multiple R-squared: 0.7201, Adjusted R-squared: 0.6423
F-statistic: 9.26 on 5 and 18 DF, p-value: 0.0001686
anova(Tol_Sym_M4_lm)
Analysis of Variance Table
Response: Sym.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.19101 0.19101 7.9874 0.011190 *
Genotype 2 0.76080 0.38040 15.9068 0.000105 ***
Site:Genotype 2 0.15542 0.07771 3.2496 0.062386 .
Residuals 18 0.43046 0.02391
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_Sym_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.12 | [0.00, 1.00]
Genotype | 0.49 | [0.18, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_Sym_M4_lm.res<-data.frame(anova(Tol_Sym_M4_lm))
Tol_Sym_M4_lm.res$Predictor<-rownames(Tol_Sym_M4_lm.res)
Tol_Sym_M4_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M4_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res<-Tol_Sym_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.863 0.0773 18 0.701 1.026
AC12 1.000 0.0773 18 0.838 1.162
AC8 0.554 0.0773 18 0.392 0.717
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.880 0.0773 18 0.717 1.042
AC12 0.622 0.0773 18 0.460 0.785
AC8 0.380 0.0773 18 0.218 0.543
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 -0.137 0.109 18 -1.249 0.4409
AC10 - AC8 0.309 0.109 18 2.829 0.0285
AC12 - AC8 0.446 0.109 18 4.078 0.0019
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.258 0.109 18 2.356 0.0734
AC10 - AC8 0.499 0.109 18 4.565 0.0007
AC12 - AC8 0.242 0.109 18 2.210 0.0965
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.863 0.0773 18 0.701 1.026
SS 0.880 0.0773 18 0.717 1.042
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 1.000 0.0773 18 0.838 1.162
SS 0.622 0.0773 18 0.460 0.785
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.554 0.0773 18 0.392 0.717
SS 0.380 0.0773 18 0.218 0.543
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0163 0.109 18 -0.149 0.8835
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS 0.3779 0.109 18 3.456 0.0028
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.1736 0.109 18 1.588 0.1297
##Save p-values
#Genotypes within Sites
Tol_Sym_M4_lm.geno<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M4_lm.geno<-Tol_Sym_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.geno$group1<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group1, sep="_")
Tol_Sym_M4_lm.geno$group2<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_Sym_M4_lm.site<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M4_lm.site<-Tol_Sym_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.site$group1<-paste(Tol_Sym_M4_lm.site$group1, Tol_Sym_M4_lm.site$Genotype, sep="_")
Tol_Sym_M4_lm.site$group2<-paste(Tol_Sym_M4_lm.site$group2, Tol_Sym_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_Sym_M4_lm.p<-rbind(Tol_Sym_M4_lm.geno[,c(1:2,4:8)], Tol_Sym_M4_lm.site[,c(1:2,4:8)])
Tol_Sym_M4_lm.p<-Tol_Sym_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_Sym_M4_lm.p$Sig<-ifelse(Tol_Sym_M4_lm.p$p<0.001, "***", ifelse(Tol_Sym_M4_lm.p$p<0.01, "**", ifelse(Tol_Sym_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_Sym_M4_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.p))
Tol_Sym_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_Sym_M4_SG<-summarySE(TolData_M4, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_Sym_M4_SG.plot<-ggplot(Tol_Sym_M4_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Symbiont Retention")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_Sym_M4_lm.p, y.position=1.1, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M4_SG.plot
##Check normality
hist(TolData$Score_Full.prop)
shapiro.test(TolData$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData$Score_Full.prop
W = 0.90488, p-value = 6.821e-05
#Not Normal
##Try square transformation
hist((TolData$Score_Full.prop)^2)
shapiro.test((TolData$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Score_Full.prop)^2
W = 0.9194, p-value = 0.0002754
#Not normal
##Try cubed transformation
hist((TolData$Score_Full.prop)^3)
shapiro.test((TolData$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Score_Full.prop)^3
W = 0.91645, p-value = 0.0002055
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreF_lm<-lm(Score_Full.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_lm)); qqline(resid(Tol_ScoreF_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_lm), resid(Tol_ScoreF_lm))
Residuals are not great, need to check for other modeling options for Color Full Set.
#Model Results
summary(Tol_ScoreF_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.217253 -0.020809 0.001419 0.037990 0.181147
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.816721 0.010646 76.716 < 2e-16 ***
Site.L 0.068538 0.015014 4.565 2.85e-05 ***
Genotype.L -0.162107 0.018189 -8.912 2.92e-12 ***
Genotype.Q -0.000730 0.018686 -0.039 0.96898
TimeP.L 0.143092 0.018189 7.867 1.44e-10 ***
TimeP.Q -0.021653 0.018686 -1.159 0.25157
Site.L:Genotype.L 0.033429 0.025724 1.300 0.19918
Site.L:Genotype.Q -0.007619 0.026324 -0.289 0.77333
Site.L:TimeP.L -0.039262 0.025724 -1.526 0.13267
Site.L:TimeP.Q -0.005353 0.026324 -0.203 0.83961
Genotype.L:TimeP.L 0.055269 0.031703 1.743 0.08686 .
Genotype.Q:TimeP.L 0.032487 0.031306 1.038 0.30393
Genotype.L:TimeP.Q -0.003858 0.031306 -0.123 0.90237
Genotype.Q:TimeP.Q -0.108144 0.033392 -3.239 0.00204 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.08798 on 55 degrees of freedom
Multiple R-squared: 0.7679, Adjusted R-squared: 0.7131
F-statistic: 14 on 13 and 55 DF, p-value: 4.691e-13
anova(Tol_ScoreF_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.18599 0.185989 24.0291 8.773e-06 ***
Genotype 2 0.58580 0.292900 37.8415 4.613e-11 ***
TimeP 2 0.48992 0.244959 31.6477 7.135e-10 ***
Site:Genotype 2 0.01497 0.007484 0.9669 0.38664
Site:TimeP 2 0.02075 0.010376 1.3405 0.27013
Genotype:TimeP 4 0.11112 0.027780 3.5891 0.01137 *
Residuals 55 0.42571 0.007740
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 0.10 | [0.01, 1.00]
Genotype | 0.32 | [0.15, 1.00]
TimeP | 0.27 | [0.10, 1.00]
Site:Genotype | 8.16e-03 | [0.00, 1.00]
Site:TimeP | 0.01 | [0.00, 1.00]
Genotype:TimeP | 0.06 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_lm.res<-data.frame(anova(Tol_ScoreF_lm))
Tol_ScoreF_lm.res$Predictor<-rownames(Tol_ScoreF_lm.res)
Tol_ScoreF_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_lm.res))
Tol_ScoreF_lm.res<-Tol_ScoreF_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Score_Full.prop)
shapiro.test(TolData_W2$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_W2$Score_Full.prop
W = 0.96098, p-value = 0.4834
#Normal
##Model as a function of Site and Genotype
Tol_ScoreF_W2_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_W2_lm)); qqline(resid(Tol_ScoreF_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_W2_lm), resid(Tol_ScoreF_W2_lm))
#Model Results
summary(Tol_ScoreF_W2_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.20185 -0.04952 -0.00045 0.04135 0.16725
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.706208 0.023330 30.270 3.15e-16 ***
Site.L 0.093421 0.032994 2.831 0.011517 *
Genotype.L -0.203806 0.040938 -4.978 0.000115 ***
Genotype.Q -0.068453 0.039874 -1.717 0.104193
Site.L:Genotype.L 0.003375 0.057895 0.058 0.954193
Site.L:Genotype.Q 0.010407 0.056391 0.185 0.855768
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1112 on 17 degrees of freedom
Multiple R-squared: 0.6882, Adjusted R-squared: 0.5964
F-statistic: 7.503 on 5 and 17 DF, p-value: 7e-04
anova(Tol_ScoreF_W2_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.12458 0.124577 10.0661 0.0055647 **
Genotype 2 0.33925 0.169625 13.7060 0.0002851 ***
Site:Genotype 2 0.00045 0.000226 0.0183 0.9819098
Residuals 17 0.21039 0.012376
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 0.18 | [0.00, 1.00]
Genotype | 0.50 | [0.17, 1.00]
Site:Genotype | 6.70e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_W2_lm.res<-data.frame(anova(Tol_ScoreF_W2_lm))
Tol_ScoreF_W2_lm.res$Predictor<-rownames(Tol_ScoreF_W2_lm.res)
Tol_ScoreF_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_W2_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res<-Tol_ScoreF_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.755 0.0556 17 0.638 0.872
AC12 0.702 0.0556 17 0.585 0.819
AC8 0.463 0.0556 17 0.346 0.581
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.890 0.0556 17 0.772 1.007
AC12 0.822 0.0556 17 0.705 0.940
AC8 0.605 0.0642 17 0.469 0.740
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0529 0.0787 17 0.673 0.7820
AC10 - AC8 0.2916 0.0787 17 3.707 0.0047
AC12 - AC8 0.2387 0.0787 17 3.034 0.0195
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0676 0.0787 17 0.859 0.6723
AC10 - AC8 0.2848 0.0850 17 3.353 0.0100
AC12 - AC8 0.2172 0.0850 17 2.557 0.0508
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.755 0.0556 17 0.638 0.872
SS 0.890 0.0556 17 0.772 1.007
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.702 0.0556 17 0.585 0.819
SS 0.822 0.0556 17 0.705 0.940
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.463 0.0556 17 0.346 0.581
SS 0.605 0.0642 17 0.469 0.740
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.135 0.0787 17 -1.713 0.1049
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.120 0.0787 17 -1.527 0.1452
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.141 0.0850 17 -1.665 0.1142
##Save p-values
#Genotypes within Sites
Tol_ScoreF_W2_lm.geno<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_W2_lm.geno<-Tol_ScoreF_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.geno$group1<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group1, sep="_")
Tol_ScoreF_W2_lm.geno$group2<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_W2_lm.site<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_W2_lm.site<-Tol_ScoreF_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.site$group1<-paste(Tol_ScoreF_W2_lm.site$group1, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
Tol_ScoreF_W2_lm.site$group2<-paste(Tol_ScoreF_W2_lm.site$group2, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_W2_lm.p<-rbind(Tol_ScoreF_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreF_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreF_W2_lm.p<-Tol_ScoreF_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_W2_lm.p$Sig<-ifelse(Tol_ScoreF_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_W2_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.p))
Tol_ScoreF_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_W2_SG<-summarySE(TolData_W2, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_W2_SG.plot<-ggplot(Tol_ScoreF_W2_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention Full Set")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_W2_lm.p, y.position=0.9, step.increase=0.35, label="Sig", hide.ns=TRUE); Tol_ScoreF_W2_SG.plot
##Check normality
hist(TolData_M1$Score_Full.prop)
shapiro.test(TolData_M1$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_M1$Score_Full.prop
W = 0.89042, p-value = 0.01921
#Not normal
##Try square transformation
hist((TolData_M1$Score_Full.prop)^2)
shapiro.test((TolData_M1$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M1$Score_Full.prop)^2
W = 0.89427, p-value = 0.02287
#Not normal
##Try cubed transformation
hist((TolData_M1$Score_Full.prop)^3)
shapiro.test((TolData_M1$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M1$Score_Full.prop)^3
W = 0.88949, p-value = 0.01842
#Not normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M1_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_M1_lm)); qqline(resid(Tol_ScoreF_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M1_lm), resid(Tol_ScoreF_M1_lm))
#Model Results
summary(Tol_ScoreF_M1_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.205567 -0.027956 0.001925 0.027094 0.192833
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.83440 0.02043 40.846 < 2e-16 ***
Site.L 0.07153 0.02889 2.476 0.024843 *
Genotype.L -0.15896 0.03357 -4.736 0.000224 ***
Genotype.Q 0.08757 0.03711 2.360 0.031325 *
Site.L:Genotype.L 0.07642 0.04747 1.610 0.126961
Site.L:Genotype.Q 0.01094 0.05248 0.208 0.837562
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.09494 on 16 degrees of freedom
Multiple R-squared: 0.6982, Adjusted R-squared: 0.6039
F-statistic: 7.403 on 5 and 16 DF, p-value: 0.0009114
anova(Tol_ScoreF_M1_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.057569 0.057569 6.3867 0.0224118 *
Genotype 2 0.252333 0.126167 13.9969 0.0003061 ***
Site:Genotype 2 0.023755 0.011877 1.3177 0.2952999
Residuals 16 0.144222 0.009014
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.12 | [0.00, 1.00]
Genotype | 0.53 | [0.19, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_M1_lm.res<-data.frame(anova(Tol_ScoreF_M1_lm))
Tol_ScoreF_M1_lm.res$Predictor<-rownames(Tol_ScoreF_M1_lm.res)
Tol_ScoreF_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M1_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res<-Tol_ScoreF_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.967 0.0475 16 0.866 1.068
AC12 0.719 0.0548 16 0.602 0.835
AC8 0.666 0.0475 16 0.565 0.766
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.998 0.0475 16 0.897 1.099
AC12 0.807 0.0548 16 0.691 0.923
AC8 0.850 0.0475 16 0.749 0.950
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.2484 0.0725 16 3.425 0.0092
AC10 - AC8 0.3012 0.0671 16 4.487 0.0010
AC12 - AC8 0.0528 0.0725 16 0.729 0.7504
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1909 0.0725 16 2.633 0.0452
AC10 - AC8 0.1484 0.0671 16 2.210 0.0998
AC12 - AC8 -0.0425 0.0725 16 -0.587 0.8292
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.967 0.0475 16 0.866 1.068
SS 0.998 0.0475 16 0.897 1.099
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.719 0.0548 16 0.602 0.835
SS 0.807 0.0548 16 0.691 0.923
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.666 0.0475 16 0.565 0.766
SS 0.850 0.0475 16 0.749 0.950
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0311 0.0671 16 -0.463 0.6499
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0885 0.0775 16 -1.142 0.2702
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1839 0.0671 16 -2.739 0.0146
##Save p-values
#Genotypes within Sites
Tol_ScoreF_M1_lm.geno<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M1_lm.geno<-Tol_ScoreF_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.geno$group1<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group1, sep="_")
Tol_ScoreF_M1_lm.geno$group2<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_M1_lm.site<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M1_lm.site<-Tol_ScoreF_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.site$group1<-paste(Tol_ScoreF_M1_lm.site$group1, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
Tol_ScoreF_M1_lm.site$group2<-paste(Tol_ScoreF_M1_lm.site$group2, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_M1_lm.p<-rbind(Tol_ScoreF_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M1_lm.p<-Tol_ScoreF_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_M1_lm.p$Sig<-ifelse(Tol_ScoreF_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_M1_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.p))
Tol_ScoreF_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_M1_SG<-summarySE(TolData_M1, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_M1_SG.plot<-ggplot(Tol_ScoreF_M1_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention Full Set")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M1_lm.p, y.position=1.1, step.increase=0.2, label="Sig", hide.ns=TRUE); Tol_ScoreF_M1_SG.plot
##Check normality
hist(TolData_M4$Score_Full.prop)
shapiro.test(TolData_M4$Score_Full.prop)
Shapiro-Wilk normality test
data: TolData_M4$Score_Full.prop
W = 0.83109, p-value = 0.0009936
#Not Normal
##Try square transformation
hist((TolData_M4$Score_Full.prop)^2)
shapiro.test((TolData_M4$Score_Full.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Score_Full.prop)^2
W = 0.8317, p-value = 0.001019
#Not Normal
##Try cubed transformation
hist((TolData_M4$Score_Full.prop)^3)
shapiro.test((TolData_M4$Score_Full.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Score_Full.prop)^3
W = 0.83161, p-value = 0.001015
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M4_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreF_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreF_M4_lm)); qqline(resid(Tol_ScoreF_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M4_lm), resid(Tol_ScoreF_M4_lm))
#Model Results
summary(Tol_ScoreF_M4_lm)
Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.137800 -0.019025 0.003675 0.020575 0.099500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.90906 0.01100 82.614 < 2e-16 ***
Site.L 0.03859 0.01556 2.480 0.0233 *
Genotype.L -0.12460 0.01906 -6.538 3.82e-06 ***
Genotype.Q -0.02191 0.01906 -1.149 0.2654
Site.L:Genotype.L 0.01901 0.02695 0.705 0.4896
Site.L:Genotype.Q -0.04168 0.02695 -1.546 0.1394
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05391 on 18 degrees of freedom
Multiple R-squared: 0.7468, Adjusted R-squared: 0.6765
F-statistic: 10.62 on 5 and 18 DF, p-value: 7.153e-05
anova(Tol_ScoreF_M4_lm)
Analysis of Variance Table
Response: Score_Full.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.017871 0.017871 6.1496 0.02326 *
Genotype 2 0.128043 0.064021 22.0308 1.452e-05 ***
Site:Genotype 2 0.008394 0.004197 1.4442 0.26199
Residuals 18 0.052308 0.002906
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.09 | [0.00, 1.00]
Genotype | 0.62 | [0.34, 1.00]
Site:Genotype | 0.04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreF_M4_lm.res<-data.frame(anova(Tol_ScoreF_M4_lm))
Tol_ScoreF_M4_lm.res$Predictor<-rownames(Tol_ScoreF_M4_lm.res)
Tol_ScoreF_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M4_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res<-Tol_ScoreF_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.982 0.027 18 0.926 1.039
AC12 0.876 0.027 18 0.819 0.932
AC8 0.787 0.027 18 0.731 0.844
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.994 0.027 18 0.937 1.051
AC12 0.978 0.027 18 0.922 1.035
AC8 0.837 0.027 18 0.780 0.893
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1069 0.0381 18 2.804 0.0301
AC10 - AC8 0.1952 0.0381 18 5.122 0.0002
AC12 - AC8 0.0883 0.0381 18 2.318 0.0788
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0157 0.0381 18 0.411 0.9115
AC10 - AC8 0.1572 0.0381 18 4.124 0.0018
AC12 - AC8 0.1415 0.0381 18 3.713 0.0043
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.982 0.027 18 0.926 1.039
SS 0.994 0.027 18 0.937 1.051
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.876 0.027 18 0.819 0.932
SS 0.978 0.027 18 0.922 1.035
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.787 0.027 18 0.731 0.844
SS 0.837 0.027 18 0.780 0.893
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0115 0.0381 18 -0.302 0.7663
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.1027 0.0381 18 -2.694 0.0148
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0495 0.0381 18 -1.299 0.2103
##Save p-values
#Genotypes within Sites
Tol_ScoreF_M4_lm.geno<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M4_lm.geno<-Tol_ScoreF_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.geno$group1<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group1, sep="_")
Tol_ScoreF_M4_lm.geno$group2<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreF_M4_lm.site<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M4_lm.site<-Tol_ScoreF_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.site$group1<-paste(Tol_ScoreF_M4_lm.site$group1, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
Tol_ScoreF_M4_lm.site$group2<-paste(Tol_ScoreF_M4_lm.site$group2, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreF_M4_lm.p<-rbind(Tol_ScoreF_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M4_lm.p<-Tol_ScoreF_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreF_M4_lm.p$Sig<-ifelse(Tol_ScoreF_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreF_M4_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.p))
Tol_ScoreF_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreF_M4_SG<-summarySE(TolData_M4, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreF_M4_SG.plot<-ggplot(Tol_ScoreF_M4_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention Full Set")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreF_M4_lm.p, y.position=1.1, step.increase=0.22, label="Sig", hide.ns=TRUE); Tol_ScoreF_M4_SG.plot
##Check normality
hist(TolData$Score_TP.prop)
shapiro.test(TolData$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData$Score_TP.prop
W = 0.9161, p-value = 0.0001986
#Not Normal
##Try square transformation
hist((TolData$Score_TP.prop)^2)
shapiro.test((TolData$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData$Score_TP.prop)^2
W = 0.92506, p-value = 0.0004895
#Not normal
##Try cubed transformation
hist((TolData$Score_TP.prop)^3)
shapiro.test((TolData$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData$Score_TP.prop)^3
W = 0.92936, p-value = 0.0007674
#Not normal
##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreTP_lm<-lm(Score_TP.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_lm)); qqline(resid(Tol_ScoreTP_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_lm), resid(Tol_ScoreTP_lm))
Residuals are not great, need to check for other modeling options for Color by Timepoint.
#Model Results
summary(Tol_ScoreTP_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + TimeP + Site:Genotype +
Site:TimeP + Genotype:TimeP, data = TolData)
Residuals:
Min 1Q Median 3Q Max
-0.133294 -0.019060 0.000725 0.026225 0.122906
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.883494 0.006319 139.809 < 2e-16 ***
Site.L 0.038546 0.008912 4.325 6.47e-05 ***
Genotype.L -0.102221 0.010797 -9.468 3.81e-13 ***
Genotype.Q 0.004971 0.011092 0.448 0.655784
TimeP.L 0.068495 0.010797 6.344 4.45e-08 ***
TimeP.Q 0.021265 0.011092 1.917 0.060411 .
Site.L:Genotype.L 0.014004 0.015269 0.917 0.363050
Site.L:Genotype.Q -0.008574 0.015625 -0.549 0.585424
Site.L:TimeP.L -0.036792 0.015269 -2.410 0.019344 *
Site.L:TimeP.Q -0.004886 0.015625 -0.313 0.755697
Genotype.L:TimeP.L 0.029281 0.018818 1.556 0.125448
Genotype.Q:TimeP.L 0.024642 0.018582 1.326 0.190292
Genotype.L:TimeP.Q 0.015340 0.018582 0.826 0.412653
Genotype.Q:TimeP.Q -0.072757 0.019821 -3.671 0.000548 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05222 on 55 degrees of freedom
Multiple R-squared: 0.763, Adjusted R-squared: 0.7069
F-statistic: 13.62 on 13 and 55 DF, p-value: 8.098e-13
anova(Tol_ScoreTP_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.057199 0.057199 20.9741 2.709e-05 ***
Genotype 2 0.238412 0.119206 43.7110 4.331e-12 ***
TimeP 2 0.117103 0.058552 21.4700 1.284e-07 ***
Site:Genotype 2 0.003221 0.001611 0.5906 0.557481
Site:TimeP 2 0.017283 0.008642 3.1687 0.049832 *
Genotype:TimeP 4 0.049580 0.012395 4.5450 0.003041 **
Residuals 55 0.149993 0.002727
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
----------------------------------------
Site | 0.09 | [0.01, 1.00]
Genotype | 0.38 | [0.20, 1.00]
TimeP | 0.19 | [0.04, 1.00]
Site:Genotype | 5.09e-03 | [0.00, 1.00]
Site:TimeP | 0.03 | [0.00, 1.00]
Genotype:TimeP | 0.08 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_lm.res<-data.frame(anova(Tol_ScoreTP_lm))
Tol_ScoreTP_lm.res$Predictor<-rownames(Tol_ScoreTP_lm.res)
Tol_ScoreTP_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_lm.res))
Tol_ScoreTP_lm.res<-Tol_ScoreTP_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Strong influence of Timepoint, including interactions with variables of interest, Site and Genotype. Will analyze each Timepoint individually.
##Check normality
hist(TolData_W2$Score_TP.prop)
shapiro.test(TolData_W2$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_W2$Score_TP.prop
W = 0.95832, p-value = 0.4302
#Normal
##Model as a function of Site and Genotype
Tol_ScoreTP_W2_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_W2)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_W2_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_W2_lm)); qqline(resid(Tol_ScoreTP_W2_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_W2_lm), resid(Tol_ScoreTP_W2_lm))
#Model Results
summary(Tol_ScoreTP_W2_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_W2)
Residuals:
Min 1Q Median 3Q Max
-0.117075 -0.028442 -0.002075 0.023658 0.099425
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.843522 0.013080 64.488 < 2e-16 ***
Site.L 0.062257 0.018498 3.366 0.00367 **
Genotype.L -0.117129 0.022952 -5.103 8.83e-05 ***
Genotype.Q -0.042426 0.022356 -1.898 0.07485 .
Site.L:Genotype.L -0.002071 0.032459 -0.064 0.94987
Site.L:Genotype.Q 0.004044 0.031616 0.128 0.89972
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.06237 on 17 degrees of freedom
Multiple R-squared: 0.7169, Adjusted R-squared: 0.6336
F-statistic: 8.61 on 5 and 17 DF, p-value: 0.0003242
anova(Tol_ScoreTP_W2_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.054151 0.054151 13.9200 0.0016623 **
Genotype 2 0.113242 0.056621 14.5549 0.0002073 ***
Site:Genotype 2 0.000083 0.000041 0.0106 0.9894534
Residuals 17 0.066133 0.003890
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
---------------------------------------
Site | 0.23 | [0.01, 1.00]
Genotype | 0.48 | [0.15, 1.00]
Site:Genotype | 3.53e-04 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_W2_lm.res<-data.frame(anova(Tol_ScoreTP_W2_lm))
Tol_ScoreTP_W2_lm.res$Predictor<-rownames(Tol_ScoreTP_W2_lm.res)
Tol_ScoreTP_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_W2_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res<-Tol_ScoreTP_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.863 0.0312 17 0.797 0.929
AC12 0.836 0.0312 17 0.771 0.902
AC8 0.699 0.0312 17 0.633 0.765
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.955 0.0312 17 0.889 1.021
AC12 0.920 0.0312 17 0.854 0.986
AC8 0.788 0.0360 17 0.712 0.864
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0263 0.0441 17 0.597 0.8236
AC10 - AC8 0.1636 0.0441 17 3.709 0.0047
AC12 - AC8 0.1373 0.0441 17 3.112 0.0166
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0354 0.0441 17 0.803 0.7065
AC10 - AC8 0.1677 0.0476 17 3.521 0.0070
AC12 - AC8 0.1323 0.0476 17 2.778 0.0328
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.863 0.0312 17 0.797 0.929
SS 0.955 0.0312 17 0.889 1.021
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.836 0.0312 17 0.771 0.902
SS 0.920 0.0312 17 0.854 0.986
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.699 0.0312 17 0.633 0.765
SS 0.788 0.0360 17 0.712 0.864
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0925 0.0441 17 -2.096 0.0513
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0834 0.0441 17 -1.890 0.0759
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.0883 0.0476 17 -1.854 0.0812
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_W2_lm.geno<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_W2_lm.geno<-Tol_ScoreTP_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.geno$group1<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group1, sep="_")
Tol_ScoreTP_W2_lm.geno$group2<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_W2_lm.site<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_W2_lm.site<-Tol_ScoreTP_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.site$group1<-paste(Tol_ScoreTP_W2_lm.site$group1, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
Tol_ScoreTP_W2_lm.site$group2<-paste(Tol_ScoreTP_W2_lm.site$group2, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_W2_lm.p<-rbind(Tol_ScoreTP_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_W2_lm.p<-Tol_ScoreTP_W2_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_W2_lm.p$Sig<-ifelse(Tol_ScoreTP_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_W2_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_W2_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.p))
Tol_ScoreTP_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_W2_SG<-summarySE(TolData_W2, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_W2_SG.plot<-ggplot(Tol_ScoreTP_W2_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention by Timepoint")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_W2_lm.p, y.position=1, step.increase=0.45, label="Sig", hide.ns=TRUE); Tol_ScoreTP_W2_SG.plot
##Check normality
hist(TolData_M1$Score_TP.prop)
shapiro.test(TolData_M1$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_M1$Score_TP.prop
W = 0.9028, p-value = 0.03383
#Not normal
##Try square transformation
hist((TolData_M1$Score_TP.prop)^2)
shapiro.test((TolData_M1$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M1$Score_TP.prop)^2
W = 0.90185, p-value = 0.03238
#Not normal
##Try cubed transformation
hist((TolData_M1$Score_TP.prop)^3)
shapiro.test((TolData_M1$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M1$Score_TP.prop)^3
W = 0.89818, p-value = 0.02735
#Not normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreTP_M1_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M1)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_M1_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M1_lm)); qqline(resid(Tol_ScoreTP_M1_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M1_lm), resid(Tol_ScoreTP_M1_lm))
#Model Results
summary(Tol_ScoreTP_M1_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M1)
Residuals:
Min 1Q Median 3Q Max
-0.125200 -0.020725 -0.001675 0.015662 0.131000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.866131 0.012681 68.303 < 2e-16 ***
Site.L 0.041582 0.017933 2.319 0.034 *
Genotype.L -0.114746 0.020837 -5.507 4.78e-05 ***
Genotype.Q 0.064377 0.023036 2.795 0.013 *
Site.L:Genotype.L 0.050375 0.029467 1.710 0.107
Site.L:Genotype.Q 0.004277 0.032577 0.131 0.897
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05893 on 16 degrees of freedom
Multiple R-squared: 0.7446, Adjusted R-squared: 0.6647
F-statistic: 9.327 on 5 and 16 DF, p-value: 0.0002606
anova(Tol_ScoreTP_M1_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.019311 0.019311 5.5598 0.03144 *
Genotype 2 0.132460 0.066230 19.0683 5.822e-05 ***
Site:Genotype 2 0.010210 0.005105 1.4698 0.25940
Residuals 16 0.055573 0.003473
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.09 | [0.00, 1.00]
Genotype | 0.61 | [0.30, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_M1_lm.res<-data.frame(anova(Tol_ScoreTP_M1_lm))
Tol_ScoreTP_M1_lm.res$Predictor<-rownames(Tol_ScoreTP_M1_lm.res)
Tol_ScoreTP_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M1_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res<-Tol_ScoreTP_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.968 0.0295 16 0.906 1.031
AC12 0.787 0.0340 16 0.715 0.859
AC8 0.755 0.0295 16 0.693 0.818
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.979 0.0295 16 0.917 1.041
AC12 0.841 0.0340 16 0.768 0.913
AC8 0.867 0.0295 16 0.805 0.930
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1815 0.0450 16 4.032 0.0026
AC10 - AC8 0.2127 0.0417 16 5.103 0.0003
AC12 - AC8 0.0312 0.0450 16 0.693 0.7710
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.1385 0.0450 16 3.077 0.0187
AC10 - AC8 0.1119 0.0417 16 2.685 0.0408
AC12 - AC8 -0.0266 0.0450 16 -0.591 0.8269
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.968 0.0295 16 0.906 1.031
SS 0.979 0.0295 16 0.917 1.041
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.787 0.0340 16 0.715 0.859
SS 0.841 0.0340 16 0.768 0.913
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.755 0.0295 16 0.693 0.818
SS 0.867 0.0295 16 0.805 0.930
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.0109 0.0417 16 -0.262 0.7970
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.0539 0.0481 16 -1.119 0.2795
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS -0.1116 0.0417 16 -2.679 0.0165
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_M1_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M1_lm.geno<-Tol_ScoreTP_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.geno$group1<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group1, sep="_")
Tol_ScoreTP_M1_lm.geno$group2<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_M1_lm.site<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M1_lm.site<-Tol_ScoreTP_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.site$group1<-paste(Tol_ScoreTP_M1_lm.site$group1, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
Tol_ScoreTP_M1_lm.site$group2<-paste(Tol_ScoreTP_M1_lm.site$group2, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_M1_lm.p<-rbind(Tol_ScoreTP_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M1_lm.p<-Tol_ScoreTP_M1_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_M1_lm.p$Sig<-ifelse(Tol_ScoreTP_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M1_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_M1_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.p))
Tol_ScoreTP_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_M1_SG<-summarySE(TolData_M1, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_M1_SG.plot<-ggplot(Tol_ScoreTP_M1_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention by Timepoint")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M1_lm.p, y.position=1.05, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M1_SG.plot
##Check normality
hist(TolData_M4$Score_TP.prop)
shapiro.test(TolData_M4$Score_TP.prop)
Shapiro-Wilk normality test
data: TolData_M4$Score_TP.prop
W = 0.85397, p-value = 0.002593
#Not Normal
##Try square transformation
hist((TolData_M4$Score_TP.prop)^2)
shapiro.test((TolData_M4$Score_TP.prop)^2)
Shapiro-Wilk normality test
data: (TolData_M4$Score_TP.prop)^2
W = 0.85439, p-value = 0.00264
#Not Normal
##Try cubed transformation
hist((TolData_M4$Score_TP.prop)^3)
shapiro.test((TolData_M4$Score_TP.prop)^3)
Shapiro-Wilk normality test
data: (TolData_M4$Score_TP.prop)^3
W = 0.85457, p-value = 0.002661
#Not Normal
##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreTP_M4_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M4)
##Check Normality of Residuals
#Distribution
plot(density(resid(Tol_ScoreTP_M4_lm)))
#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M4_lm)); qqline(resid(Tol_ScoreTP_M4_lm))
##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M4_lm), resid(Tol_ScoreTP_M4_lm))
#Model Results
summary(Tol_ScoreTP_M4_lm)
Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype,
data = TolData_M4)
Residuals:
Min 1Q Median 3Q Max
-0.070400 -0.009438 0.004100 0.011281 0.057100
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.940608 0.006243 150.654 < 2e-16 ***
Site.L 0.010536 0.008830 1.193 0.2483
Genotype.L -0.075254 0.010814 -6.959 1.68e-06 ***
Genotype.Q -0.007308 0.010814 -0.676 0.5078
Site.L:Genotype.L -0.006950 0.015293 -0.454 0.6549
Site.L:Genotype.Q -0.032086 0.015293 -2.098 0.0503 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.03059 on 18 degrees of freedom
Multiple R-squared: 0.7531, Adjusted R-squared: 0.6846
F-statistic: 10.98 on 5 and 18 DF, p-value: 5.763e-05
anova(Tol_ScoreTP_M4_lm)
Analysis of Variance Table
Response: Score_TP.prop
Df Sum Sq Mean Sq F value Pr(>F)
Site 1 0.001332 0.0013321 1.4238 0.2483
Genotype 2 0.045732 0.0228662 24.4414 7.407e-06 ***
Site:Genotype 2 0.004311 0.0021557 2.3042 0.1285
Residuals 18 0.016840 0.0009355
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)
# Effect Size for ANOVA (Type I)
Parameter | Eta2 | 95% CI
-----------------------------------
Site | 0.02 | [0.00, 1.00]
Genotype | 0.67 | [0.41, 1.00]
Site:Genotype | 0.06 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].
##Save model results
Tol_ScoreTP_M4_lm.res<-data.frame(anova(Tol_ScoreTP_M4_lm))
Tol_ScoreTP_M4_lm.res$Predictor<-rownames(Tol_ScoreTP_M4_lm.res)
Tol_ScoreTP_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M4_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res<-Tol_ScoreTP_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")
Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.
#Pairwise comparisons across:
#Genotypes within Sites
emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)
$emmeans
Site = KL:
Genotype emmean SE df lower.CL upper.CL
AC10 0.989 0.0153 18 0.957 1.021
AC12 0.921 0.0153 18 0.888 0.953
AC8 0.890 0.0153 18 0.858 0.922
Site = SS:
Genotype emmean SE df lower.CL upper.CL
AC10 0.993 0.0153 18 0.960 1.025
AC12 0.973 0.0153 18 0.940 1.005
AC8 0.879 0.0153 18 0.847 0.911
Confidence level used: 0.95
$contrasts
Site = KL:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0686 0.0216 18 3.171 0.0139
AC10 - AC8 0.0995 0.0216 18 4.599 0.0006
AC12 - AC8 0.0309 0.0216 18 1.429 0.3478
Site = SS:
contrast estimate SE df t.ratio p.value
AC10 - AC12 0.0199 0.0216 18 0.922 0.6335
AC10 - AC8 0.1134 0.0216 18 5.242 0.0002
AC12 - AC8 0.0934 0.0216 18 4.320 0.0011
P value adjustment: tukey method for comparing a family of 3 estimates
#Sites within Genotypes
emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)
$emmeans
Genotype = AC10:
Site emmean SE df lower.CL upper.CL
KL 0.989 0.0153 18 0.957 1.021
SS 0.993 0.0153 18 0.960 1.025
Genotype = AC12:
Site emmean SE df lower.CL upper.CL
KL 0.921 0.0153 18 0.888 0.953
SS 0.973 0.0153 18 0.940 1.005
Genotype = AC8:
Site emmean SE df lower.CL upper.CL
KL 0.890 0.0153 18 0.858 0.922
SS 0.879 0.0153 18 0.847 0.911
Confidence level used: 0.95
$contrasts
Genotype = AC10:
contrast estimate SE df t.ratio p.value
KL - SS -0.00332 0.0216 18 -0.154 0.8795
Genotype = AC12:
contrast estimate SE df t.ratio p.value
KL - SS -0.05195 0.0216 18 -2.402 0.0273
Genotype = AC8:
contrast estimate SE df t.ratio p.value
KL - SS 0.01057 0.0216 18 0.489 0.6308
##Save p-values
#Genotypes within Sites
Tol_ScoreTP_M4_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M4_lm.geno<-Tol_ScoreTP_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.geno$group1<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group1, sep="_")
Tol_ScoreTP_M4_lm.geno$group2<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group2, sep="_")
#Sites within Genotypes
Tol_ScoreTP_M4_lm.site<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M4_lm.site<-Tol_ScoreTP_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.site$group1<-paste(Tol_ScoreTP_M4_lm.site$group1, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
Tol_ScoreTP_M4_lm.site$group2<-paste(Tol_ScoreTP_M4_lm.site$group2, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
#Full list of p-values
Tol_ScoreTP_M4_lm.p<-rbind(Tol_ScoreTP_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M4_lm.p<-Tol_ScoreTP_M4_lm.p %>% dplyr::rename( p = p.value)
#Add Significance Levels
Tol_ScoreTP_M4_lm.p$Sig<-ifelse(Tol_ScoreTP_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M4_lm.p$p<0.05, "*", NA)))
#Specify Response and Timepoint
Tol_ScoreTP_M4_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.p))
Tol_ScoreTP_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.p))
##Summary statistics by Site and Genotype
Tol_ScoreTP_M4_SG<-summarySE(TolData_M4, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)
##Plot Average Retention across Treatments
Tol_ScoreTP_M4_SG.plot<-ggplot(Tol_ScoreTP_M4_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) +
geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
geom_point(size=point.sz)+
ggtitle("Color Retention by Timepoint")+
theme_classic()+
theme( axis.title.x = element_text(size = axis.title.sz),
axis.title.y = element_text(size = axis.title.sz),
axis.text.x=element_text(size=axis.txt.sz, colour="black"),
axis.text.y=element_text(size=axis.txt.sz, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
legend.box.background = element_rect(color = "black"),
legend.position="bottom",
legend.direction="horizontal",
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="Proportion Retained")+
ylim(0, 1.5)+
scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
scale_color_manual(values = Geno.colors.o)+
stat_pvalue_manual(data=Tol_ScoreTP_M4_lm.p, y.position=1.1, step.increase=0.4, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M4_SG.plot
Creating a heatmap to compare the direction and significance of pairwise comparisons across Tolerance metrics. Positive estimates indicate that Group 1 > Group 2. P values of less than 0.05 are considered significant.
##Combine Results
Tol_contrasts<-rbind(Tol_Chl_W2_lm.p, Tol_Chl_M1_lm.p, Tol_Chl_M4_lm.p,
Tol_Sym_W2_lm.p, Tol_Sym_M1_lm.p, Tol_Sym_M4_lm.p,
Tol_ScoreF_W2_lm.p, Tol_ScoreF_M1_lm.p, Tol_ScoreF_M4_lm.p,
Tol_ScoreTP_W2_lm.p, Tol_ScoreTP_M1_lm.p, Tol_ScoreTP_M4_lm.p)
##Create Contrast Column
Tol_contrasts$Contrast<-paste(Tol_contrasts$group1, Tol_contrasts$group2, sep=" v. ")
##Add Set column with Contrast and Timepoint
Tol_contrasts$Set<-paste(Tol_contrasts$TimeP, Tol_contrasts$Contrast, sep=" ")
##Re-order by Seasonal Timepoint
Tol_contrasts$TimeP<-factor(Tol_contrasts$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
Tol_contrasts<- Tol_contrasts %>% arrange(desc(as.numeric(TimeP)))
Tol_contrasts$Set<-factor(Tol_contrasts$Set, levels=c(unique(Tol_contrasts$Set)), ordered=TRUE)
Convert Estimate to the Response scale (instead of log +1 scale) for models where the response was log transformed
Tol_contrasts$Estimate<-Tol_contrasts$estimate
#Chl M4
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")])-1)
#Sym M1
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")])-1)
Tol_Pairs_W2.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="W2"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Summer 1")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_W2.plot
Tol_Pairs_M1.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M1"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Summer 2")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_M1.plot
Tol_Pairs_M4.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M4"),], aes(x=Response, y=Contrast, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts Month 4")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_M4.plot
Tol_Pairs.plot<-ggplot(data=Tol_contrasts, aes(x=Response, y=Set, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Pairwise Contrasts")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs.plot
Subsetting Contrasts to only include contrasts with significant differences in at least one of the thermal tolerance metrics (Chlorophyll or Symbiont or Color retention)
##Remove rows with non-significant p-values
Tol_contrasts_sig<-subset(Tol_contrasts, p < 0.05)
##Keep Sets where at least one metric shows a significant difference
Tol_Sets_sig<-data.frame(Set=c(unique(Tol_contrasts_sig$Set)))
Tol_contrasts_sig_Set<-merge(Tol_Sets_sig, Tol_contrasts, all.x=TRUE, all.y=FALSE)
Tol_Pairs_sig.plot<-ggplot(data=Tol_contrasts_sig_Set, aes(x=Response, y=Set, fill=Estimate))+
geom_tile()+
geom_text(aes(label=Sig), size=sig.sz)+
scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
theme_bw()+
ggtitle("Significant Pairwise Contrasts")+
theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
legend.text=element_text(size=leg.txt.sz),
legend.title=element_text(size=leg.title.sz),
plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
labs(x="", y="", fill="Estimate");Tol_Pairs_sig.plot
##Create Panel
Tolerance_fig<-plot_grid(Tol_Chl_W2_SG.plot, Tol_Chl_M1_SG.plot, Tol_Chl_M4_SG.plot,
Tol_ScoreF_W2_SG.plot, Tol_ScoreF_M1_SG.plot, Tol_ScoreF_M4_SG.plot,
Tol_ScoreTP_W2_SG.plot, Tol_ScoreTP_M1_SG.plot, Tol_ScoreTP_M4_SG.plot,
Tol_Sym_W2_SG.plot, Tol_Sym_M1_SG.plot, Tol_Sym_M4_SG.plot,
nrow=4, ncol=3,
rel_widths=1,
rel_heights=1,
labels=c("A", "B", "C",
"D", "E", "F",
"G", "H", "I",
"J", "K", "L"),
label_size=panel.lab.sz,
label_fontface = "bold")
##Save Figure
ggsave(filename="Figures/FigureS2_Tolerance_Across_Metrics.png", plot=Tolerance_fig, dpi=300, width=18, height=24, units="in")
##Save Figure
ggsave(filename="Figures/Figure3_Tolerance_Heatmap.png", plot=Tol_Pairs.plot, dpi=300, width=8, height=12, units="in")
##Combine Results Tables
TableS2_Tol.lm.res<-rbind(Tol_Chl_W2_lm.res, Tol_Chl_M1_lm.res, Tol_Chl_M4_lm.res,
Tol_Sym_W2_lm.res, Tol_Sym_M1_lm.res, Tol_Sym_M4_lm.res,
Tol_ScoreF_W2_lm.res, Tol_ScoreF_M1_lm.res, Tol_ScoreF_M4_lm.res,
Tol_ScoreTP_W2_lm.res, Tol_ScoreTP_M1_lm.res, Tol_ScoreTP_M4_lm.res)
##Organize
names(TableS2_Tol.lm.res)
[1] "DF" "Sum.Sq" "Mean.Sq" "F.value" "p.value" "Predictor"
[7] "EtaSq" "Response" "TimeP"
TableS2_Tol.lm.res<-TableS2_Tol.lm.res[,c("TimeP", "Response", "Predictor", "DF", "Sum.Sq", "Mean.Sq", "F.value", "p.value", "EtaSq")]
#Round to 4 digits
TableS2_Tol.lm.res$Sum.Sq<-round(TableS2_Tol.lm.res$Sum.Sq, 4)
TableS2_Tol.lm.res$Mean.Sq<-round(TableS2_Tol.lm.res$Mean.Sq, 4)
TableS2_Tol.lm.res$F.value<-round(TableS2_Tol.lm.res$F.value, 4)
TableS2_Tol.lm.res$p.value<-round(TableS2_Tol.lm.res$p.value, 4)
TableS2_Tol.lm.res$EtaSq<-round(TableS2_Tol.lm.res$EtaSq, 4)
##Write Out Table
write.csv(TableS2_Tol.lm.res, "Tables/TableS2_Tolerance_LM_Results.csv", row.names=FALSE)
##Pairwise Results Table
TableS3_Tol.pairwise<-Tol_contrasts
##Organize
names(TableS3_Tol.pairwise)
[1] "group1" "group2" "estimate" "SE" "df" "t.ratio" "p"
[8] "Sig" "Response" "TimeP" "Contrast" "Set" "Estimate"
TableS3_Tol.pairwise<-TableS3_Tol.pairwise[,c("TimeP", "Response", "Contrast", "Estimate", "SE", "df", "t.ratio", "p")]
TableS3_Tol.pairwise<-TableS3_Tol.pairwise %>% dplyr::rename( DF = df) %>% dplyr::rename( p.value = p)
#Round to 4 digits
TableS3_Tol.pairwise$Estimate<-round(TableS3_Tol.pairwise$Estimate, 4)
TableS3_Tol.pairwise$SE<-round(TableS3_Tol.pairwise$SE, 4)
TableS3_Tol.pairwise$t.ratio<-round(TableS3_Tol.pairwise$t.ratio, 4)
TableS3_Tol.pairwise$p.value<-round(TableS3_Tol.pairwise$p.value, 4)
##Write Out Table
write.csv(TableS3_Tol.pairwise, "Tables/TableS3_Tolerance_Pairwise_Results.csv", row.names=FALSE)